#Tables and Figures: Quick. But Impactful? United Nations Quick Impact Projects and Violence against Civilians in Civil War
rm(list=ls())
setwd("YOUR WORKING DIRECTORY HERE")

library(tidyverse)
library(sf)
library(stargazer)
library(jtools)
library(marginaleffects)
library(ggpubr)

#FIGURE 1: Number of QIPs and Attacks against Civilians
#open shape files to create maps
car <- st_read("caf_admbnda_adm1_200k_sigcaf_reach_itos_v2/caf_admbnda_adm1_200k_sigcaf_reach_itos_v2.shp")
car$admin <- car$admin1Name
car$country <- "Central African Republic"
car <- dplyr::select(car, admin, country)

mali <- st_read("mali_admin_shp/mli_admbnda_adm2_1m_gov_20211220.shp")
mali$admin <- mali$ADM2_FR
mali$country <- "Mali"
mali <- dplyr::select(mali, admin, country)

drc <- st_read("cod_admbnda_rgc_itos_20190911_shp/cod_admbnda_adm2_rgc_20190911.shp")
drc$admin <- drc$ADM2_FR
drc$country <- "Democratic Republic of Congo"
drc <- dplyr::select(drc, admin, country)
drc <- drc %>%
  group_by(admin) %>%
  mutate(row_number = ifelse(duplicated(admin), paste0("_", row_number()), "")) %>%
  ungroup() %>%
  mutate(admin = paste0(admin, row_number)) %>%
  dplyr::select(-row_number)

ss <- st_read("ssd_admbnda_adm2_imwg_nbs_20180817/ssd_admbnda_adm2_imwg_nbs_20180817.shp")
ss$admin <- ss$ADM2_EN
ss$country <- "South Sudan"
ss <- dplyr::select(ss, admin, country)

polygons <- rbind(mali, drc, ss, car)

#data on violence against civilians
acled <- read.csv("vac_acled.csv") #subset original acled data, event_type=="Violence against civilians" & inter1==2
rebel.violence.against.civilians <- st_as_sf(acled, coords = c("longitude", "latitude"), crs = 4326) 
rebel.vac.formap <- rebel.violence.against.civilians

#open data on QIPs
df <- read.csv("DPO-QIP.csv")
#drop observations missing lon and lat
df <- subset(df, !is.na(longitude1) & !is.na(latitude1))
#drop observations with dates that are missing/do not make sense
df <- subset(df, start_date!="")
df <- subset(df, qip_ref!="2023-24 /6")
df <- subset(df, end_date!="FO to conduct final visit")
df <- subset(df, end_date!="") #drop projects with no end dates or when the end date is before the start date (which must be an error but makes project unuseable)
#drop obs where enddate is before startdate
df <- df %>%
  mutate(startdate = as.Date(start_date, format = "%Y-%m-%d"),
         enddate = as.Date(end_date, format = "%Y-%m-%d")) %>%
  filter(enddate >= startdate)
#create year and month variables
df$year <- format(as.Date(df$start_date, format="%Y-%m-%d"),"%Y")
df$month <- format(as.Date(df$start_date, format="%Y-%m-%d"),"%m")

#drop obs that have an obvious error in location
df$latitude1 <- as.numeric(df$latitude1)
df <- subset(df, latitude1<36)

#transform into spatial object
qip_start <- st_as_sf(df, coords = c("longitude1", "latitude1"), crs = 4326) 

#match QIPs to admin based on lat and lon
joined_qip.s <- st_join(qip_start, polygons, join = st_within)
qipstart <- joined_qip.s %>% 
  group_by(month, year, admin) %>%
  summarise(qipstart = n()) %>%
  ungroup()

formap <- qipstart %>%
  group_by(admin)%>%
  tally(name = "qips")
formap <- as.data.frame(formap)
map <- left_join(polygons, formap, by="admin")
rebel.vac.formap <- subset(rebel.vac.formap, year>2016 & year<2023)

create_map <- function(country_name) {
  polygons <- subset(map, country == country_name)
  points <- subset(rebel.vac.formap, country == country_name)
  
  ggplot() +
    geom_sf(data = polygons, aes(fill = qips), color = "white") +
    geom_sf(data = points, colour = 'black', size = 1, shape = 4) +
    scale_fill_gradient(low = "#0000FF33", high = "#0000FF", name = "QIPs") +
    labs(title = country_name) +
    theme_minimal()
}

car <- create_map("Central African Republic")
mali <- create_map("Mali")
drc <- create_map("Democratic Republic of Congo")
ss <- create_map("South Sudan")

ggarrange(car, mali, drc, ss)

#open data for analysis 
merge <- read.csv("qip_civilians_replication_tables_figures.csv")

#subset to only observations that had at least one attack on civilians & at least one project 
temp <- merge %>% 
  group_by(admin)%>%
  filter(sd(rebel.violence.against.civilians.lead1, na.rm = TRUE) > 0 & sd(log.cost_qipongoing, na.rm = TRUE) > 0) %>% 
  ungroup()

#TABLE 1: Descriptive statistics 
descriptives <- dplyr::select(temp, rebel.violence.against.civilians, log.qipongoing, log.cost_qipongoing, battles, gov.violence.against.civilians, mission.duration, pktroops, n_actors)
stargazer(as.data.frame(descriptives), type="html", digits=2, out="descriptives.doc")

#FIGURE 2: The Effect of QIPs on Rebel Violence against Civilians 
summary(m1 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.cost_qipongoing + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))
summary(m2 <- glm.nb(rebel.violence.against.civilians.lead1 ~ log.qipongoing + rebel.violence.against.civilians + gov.violence.against.civilians + battles + mission.duration + log.pktroops + n_actors + factor(admin), data=temp))

#jtools
plot_summs(m1, m2,
           coefs = c("QIPs-cost (logged)" = "log.cost_qipongoing", "QIPs-number (logged)" = "log.qipongoing", "Rebel violence against civilians" = "rebel.violence.against.civilians", "Government violence against civilians" = "gov.violence.against.civilians", "Battles" = "battles", "Mission duration" = "mission.duration", "PK troops (logged)" = "log.pktroops", "Number of armed actors" = "n_actors", "Other aid projects (logged)" = "log.other.aid.projects"))

#FIGURE 3: District-Month Characteristics of QIPs
summary(m1 <- glm.nb(qipongoing ~ battles, data=merge))
summary(m2 <- glm.nb(qipongoing ~ log.pktroops, data=merge))
summary(m3 <- glm.nb(qipongoing ~ pkbase, data=merge))

p1 <- plot_predictions(m1, "battles")+
  labs(y="Predicted QIPs (logged)", x="Battles", title="Violence")+
  coord_cartesian(xlim=c(0,20), ylim = c(0, 32))
p2 <- plot_predictions(m2, "log.pktroops")+
  labs(y="Predicted QIPs (logged)", x="PK troops (logged)", title="PK troops")+
  coord_cartesian(ylim = c(0, 32))
p3 <- plot_predictions(m3, "pkbase")+
  labs(y="Predicted QIPs (logged)", x="Bases", title="PK bases")+
  coord_cartesian(xlim=c(0,3),ylim = c(0, 32))

ggarrange(p1, p2, p3, nrow=1)

#TABLE 2: Percentage of District-Months With(out) QIPs and Peacekeepers (PKs)
temp$pktroops.binary.c <- "pktroops:no"
temp$pktroops.binary.c[temp$pktroops.binary==1] <- "pktroops:yes"
temp$qipongoing.binary.c <- "qip:no"
temp$qipongoing.binary.c[temp$qipongoing.binary==1] <- "qip:yes"
prop.table(table(temp$pktroops.binary.c, temp$qipongoing.binary.c))

